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DYNAMIC TRANSITIONS OF QUASI-GEOSTROPHIC CHANNEL FLOW 


HENK DIJKSTRA, TAYLAN SENGUL, JIE SEEN, AND SHOUHONG WANG 

Abstract. The main aim of this paper is to describe the dynamic transitions in flows 
described by the two-dimensional, barotropic vorticity equation in a periodic zonal channel. 
In [^, the existence of a Hopf bifurcation in this model as the Reynolds number crosses a 
critical value was proven. In this paper, we extend the results in [ 3 ] by addressing the stability 
problem of the bifurcated periodic solutions. Our main result is the explicit expression of a 
non-dimensional number 7 which controls the transition behavior. We prove that depending 
on 7, the modeled flow exhibits either a continuous (Type I) or catastrophic (Type II) 
transition. Numerical evaluation of 7 for a physically realistic region of parameter space 
suggest that a catastrophic transition is preferred in this flow. 


1 . Introduction 

Climate variability exhibits recurrent large-scale patterns which are directly linked to dy¬ 
namical processes represented in the governing dissipative dynamical system [6l [TJ [8] . The 
study of the persistence of these patterns and the transitions between them also play a crucial 
role in understanding climate change and in interpreting future climate projections [^. 

Current climate models used for developing such projections are based on the conserva¬ 
tion laws of fluid mechanics and consist of systems of nonlinear partial differential equations 
(PDEs). These can be put into the perspective of infinite-dimensional dissipative systems ex¬ 
hibiting large-dimensional attractors. The global attractor is a mathematical object strongly 
connected to the overall dissipation in the system. Climate variability is, however, often 
associated with dynamic transitions between different regimes, each represented by local 
attractors. 

There are many examples of climate phenomena where such transitions have been in¬ 
vestigated numerically, such as the transition to oscillatory behavior in models of the El 
Niho/Southern Oscillation phenomenon in the equatorial Pacific, the transitions between dif¬ 
ferent mean flow patterns of the Kuroshio Current in the North Pacific and the transitions 
between blocked and zonal flows in the midlatitude atmosphere (see, e.g., 0). However, 
rigorous mathematical results on the type of the transitions in these systems of PDEs are 
extremely scarce. 

This paper arises out of a research program to generate rigorous mathematical results on 
climate variability developed from the viewpoint of dynamical transitions. We have shown [TO] 
that the transitions of all dissipative systems can be classified into three classes: continuous, 
catastrophic and random, which correspond to very different dynamical transition behavior 
of the system. 


Key words and phrases, quasi-geostrophic flow, channel flow, spatial-temporal patterns, dynamic transi¬ 
tions, climate variability. 
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We here focus on the dynamic transitions in flows described by one of the cornerstone 
dynamical models in both oceanic and atmospheric dynamics, the two-dimensional, dimen¬ 
sionless barotropic vorticity equation given by 


( 1 . 1 ) 


T/ / A / \ rn A 2 / 

+ eJ{ iIj, Alp) + ^ = EA^ip -)- sm iry, 


where A is the Laplacian operator, J{F,G) = {dF/dx){dG/dy) — {dF/dy){dG/dx) is the 
advection operator and V' tlie geostrophic stream function. The equation (1.1) describes 


flows with a typical length scale L on a mid-latitude beta-plane with a central latitude 6q 
and a planetary vorticity gradient /3o. It can be derived from the primitive equations by the 
so-called quasi-geostrophic (QG) approximation, which assumes a dominant balance between 
the Coriolis force and the pressure gradient force m- 

To derive (O), time, length and stream function were non-dimensionalised with l/(/3oT), 
L and UL, where U is a characteristic horizontal velocity. The positive parameters e and E 
are Rossby and Ekman numbers, respectively, given by 

U ^ Ah 


e = 


E = 


where Ah is the lateral friction coefficient. The Reynolds number R is defined as 

UL 


( 1 . 2 ) 


R= ^ = 


E Ah 


The forcing term Ur sin Try in (1.1) may represent the transfer of angular momentum into 


midlatitudes due to tropical Hadley cell in an atmospheric model. In this case, the magnitude 
of the velocity can be scaled such that a,- = 1. In the ocean case, the forcing term represents 
the dimensionless wind stress 

/ p,\ 

r = —(cos Try, (J) 

TT 

Such a wind stress mimics the annually averaged zonal wind distribution over the North 
Atlantic and North Pacific with westerly (i.e. eastward) winds over the midlatitudes and 
easterlies in the tropics and polar latitudes. When the dimensionless wind stress has a 
magnitude tq and the ocean basin a depth D and the water a density p, the factor Ur is given 
by 

= '^0 

“ pDL/3oU 

In this case, we can choose U such that Or = 1, which is usually referred to as the Sverdrup 


scaling. In both ocean and atmosphere cases the equation (1.1) has two free parameters (out 
of the three e, E and R) which we choose here as R and E. 

We consider flows in a so-called zonal channel of length 2/a with walls bounding the 


domain at y = ±1 and periodic conditions in zonal direction. The equation (1.1) is therefore 
supplemented with boundary conditions 

U=0 |a;=2/a) 

/I _ 

W \y=±l Qy2 ly=±l b. 

The equations (1.1) with (1.3) admit the following steady state 

1 


(1.4) 


V'o = - 


tt^E 


sm Try, 
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which represents a zonal jet with zonal velocity field uq = —dij^o/dy = 1/{'k^E) cosvry. 

It is shown in [3] that for any a > \/3/2, V’o is both linearly and nonlinearly stable. Also, 
there is an uq with -v/3/4 < oq < V3/2 such that for any -v/3/4 < a < uq, there is a critical 
Reynolds number Rq > Q depending on a such that a simple pair of complex conjugate 
eigenvalues cross the imaginary axis as the Reynolds number R crosses Rq, leading to the 
existence of Hopf bifurcation at the critical Reynolds number. However, the stability of the 
bifurcated periodic solutions and the dynamic transition behavior near Rq are so far unknown. 
The main difficulty is caused by the lack of explicit analytical form of the eigenfunctions. 

The main objective of this article is to investigate the dynamic transition and the stability 
of the basic state (1.4) as the Reynolds number crosses a critical threshold Rq. The main result 


we obtain is that the dynamic transition from this state to new states is either continuous 
(Type-I) or catastrophic (Type-II), and is determined by the sign of a computable parameter 7 
given by (2.10). Our numerical investigations indicate that in a physically relevant parameter 


regime, only catastrophic transitions occur. 


2. Main Theorem 

Throughout Rez, Imz, z will denote the real part, imaginary part and conjugate of a 
complex number 2 . D = -^ is the derivative operator, H = (0, 2/a) x (—1,1) C and (•, •) 
is the L^(H) inner product. 

2.1. Functional setting. Considering the deviation i/' = i/ —from the basic steady state 


(1.4) and omitting the primes, we obtain from (1.1) 


( 2 . 1 ) 


T/ / A / A 2^V’x dip rn A 2 / 

at +*A</-.AV') = -;3Cos,t!,(A-+tr -)--+EAi,. 


Note that (1.1) can also be formulated in velocity u = (—§ 7 , iy) and pressure p as 


( 2 . 2 ) 


dy ’’ dx ' 

d\i 

— + e(u • V)u + /k X u = FAu — Vp + r. 


V • u = 0. 

Here / is the dimensionless Coriolis parameter on a midlatitude beta plane that gives rise to 


the dip/dx term in ( 1 . 1 ); k is the unit vector in the z-direction. 

rite the proble 

= Lu + G{u), 


(2.3) 


Using the formulation (2.2), we can write the problem in the following abstract form 

du 

dt 


where L : Hi H is the linear operator, G : Hi ^ H is the nonlinear operator, and 


H = {u = {u,v) G (L 2 (fl))^ 

Hi = HniH\n)f. 


y=±l= 0, divu = 0,U |a,=0= U U=2/a} 


The eigenvalue problem for the linearized equation of (2.1) reads 
(2.4) 


R 


EA^ip + ^ cos( 7 ry)(Ai/ 2 , + tt ipx) - ipx = IdAip 


TT"’ 


with boundary conditions (1.3). 
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Since the solution -0 is periodic in x with period 2/a, we can expand tp in Fourier series, 
so for the m-mode of the expansion, we can write 


(2.5) 


?/= am = armT. 


where Y{y) satisfies the boundary condition ^(±1) = D'^Y{±1) = 0. Plugging (2.5) into 
(2.4), we obtain a sequence of one-dimensional problem: 

R 


( 2 . 6 ) 


E{D^ - af^YY + ia 


TT^ 


cos(7r?/)(F» -a^ + TT )-i y = (3{D - a^)Y. 


The eigenvectors of (2.4) are Ym,j = E°‘""^Ym,j{y) where Ymj are the eigenvectors of (2.6) 
corresponding to the eigenvalues Pmj £ C where m G Z, j = 1, 2,3,.... Moreover f3rn,j are 
ordered so that Re f3m,j < Re (3m,k if j > k. Also (3-m,j = I3m,j and we can take Y-m,j = Ym,j- 
In particular Yo,j and (3o,j are real. 

We also need to consider the eigenvalue problem for the adjoint linear operator which can 
be written as 


R 


R 


EA Y* -7A(cos(7ry)V’*)-cos(7ry)V'* + C = l^Aip 


7^0 


TT^ 


Using "0* = we obtain the analog of (2.6) 


(2.7) 


E{D^ - a^YY* - ia. 


R 


R 


TT"’ 


D^iY* cos{'Ky)) -\—^ cos(7ry)(7r"^ — a^)Y* — Y^ 


TT"' 


= P*{D^-al)Y\ 


By basic properties of the adjoint linear eigenvalue problem, we have = f3m,j and 
(2.8) (umj,u;fc) = 0, if (m, j) / {n,k), 

where Um,j = (|;, 


2.2. The main theorem and its proof. Our main aim is to identify the transitions of 
(2.1) in the case where two complex conjugate eigenvalues cross the imaginary axis. Thus, 
we assume the following condition on the spectrum of the linearized operator. 


Assumption 1. Depending on a and E, there exists a critical Reynolds number Rq and a 
zonal wave integer mg > 1 such that 


Re(/3mo,i(-R)) = Re(/3 (R)) 


<0 if R < Ro, 
< = 0 if R = Ro, 
>0 if R > Rq- 


Re{(3m,j{Ro)) < 0, if {m,j) / (±mo, 1). 


In [3|, the validity of the Assumption is shown with mg = 1 when \/3/4 < a < og for 
some \/3/4: < ag < \/3/2 . 
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Let US define 


(2.9) 


-^1 = y Ymo,i{y)dy, 

-1 

1 

h{k) = j cofi{k'ny)Y^^^^{y) {{aniQ-Kf - - L>^) Yraa,i{y)dy, 


-1 

1 


h{k) = / sYi{k'Ky)Ymo,i{y)k^^rno,i{y)dy. 


-1 


where Ym^,! and solutions of (2.6) and (2.7) respectively for m = mo, at i? = i?o- 

We define the transition number 


( 2 . 10 ) 


7 = 


E{RQamof‘'K Im{Ii/2(fc)} 


2|/i|' 


k=l 


k 


As the next theorem shows, the sign of 7 determines the type of transition of the system 
at the critical Reynolds number Rq. 


Theorem 1. Let ^ he defined by (2.10) and let 
( 2 . 11 ) 


ifit,x,y) = + o(\/l Re/3mo,i|)- 


Ubi 


Under the Assumption^^ the following assertions hold true: 

(1) If ^ < 0 then the problem undergoes a Type-I (continuous) transition at R = Rq and 
bifurcates to the time periodic solution u^if on R > Rq which is an attractor. 

(2) // 7 > 0 then the problem undergoes a Type-II (catastrophic) transition at R = Rq 
and bifurcates to the time periodic solution u^ij on R < Rq which is a repeller. 

Proof Let u^o,! = (^,-^)i/’mo,i where V'mo,! = e“'"o’"^ymo,i(y), denote the first critical 
eigenfunction corresponding to the eigenvalue fimo,! iii Assumption 0- For simplicity of 
notation, we will denote 


= Re u 


mo,l) 


= Im u. 


mo,It 


( 2 . 12 ) 


= ReV’mo,!, = Irni/^mo,!, 


Y^ =ReYr 


mo,l) 


Y^ = ImYr 


mo,l) 


where Tmo,i solves equation (2.6). 

In the proof we will use the following trilinear operators 


(2.13) 


G{ui, uj, u;^) = -e (u/ • V)uj • UKdxdy. 
JQ 


and 


G^(u/,uj,ui^) = G(u/,uj,ui^) + G(uj,u/,uj^). 
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Step 1. Computation of nonlinear interactions. It is easy to see that only (0,fc) and {2m, k) 
adjoint modes interact nonlinearly with the critical (mo, l)-mode i.e. 

(2.14) G(u*,u^u* ;,) = 0 if n / 0, orn / 2m, = 1,2. 

We will first investigate these interactions. We have ugfe = Ug^ = (g|V’OA:,0) — (-DIofc)O) 
where Yofc is given by (2.42). In particular Ug^ is real. Using 

/■2/a 

G(u, u, u*fc) = C / = 0, 

Jx=0 

where C = fy=_i (only y-dependent terms)^//, we get 

(2.15) 0 = G(u\u\u;5fc) - Giu‘^,u^,u*Qi,)+i{G{u^,u^,uli,) + G{u‘^,u\ulk)) 

Also noting, 

, /■2/a 

G(u, u, fc) = G / e-*“-o*e*“-o^T(a2mx)da: = 0, 

Jx=0 

where G = /J=_i (only y-dependent terms)(iy and T = sin or T = cos, we get 

(2.16) 0 = G(u\ fc) + G(u2, u^, fc) -p i(G(u\ fc) - G(u^ u\ ;.)) 

Let us define 

(2.17) 


%i = G{u\u^,ulf^), 

92m,k = G{u\u\y^lm,k)- 


(2.15),(2.16) and (2.17) imply that 


(2.18) 


9ok 9f)k ; 


9f)k doki 


11 _ 22 12 _ 21 
92m^k 92m,k") 92m,k 92m,k' 


A lengthy but straightforward calculation which can also be verified by a symbolic computa¬ 
tion software shows that 

G(u\u\u^)„^,) = G(u\u2,u;^fc), 

G(u\u\u^^^,) = -G(u\u2,u^)„^,), 

which implies that 

(2.19) = igll,^k. 

Step 2. Approximation of the center manifold. Next we obtain an approximation for the 
center manifold function 4>. Let H = Ei ® E 2 , Ei = span{u^, }, E 2 = Ei. By m , near 
R = Ro^ the center manifold $ can be approximated by the formula 

{{-Cf + 4Im(/3)2)(-£)$ = ((-£)2 + 2lm{f3f)P2G{xiu^ + xsu^) 

+ 2Im(/3)^P2G(a:iu^ - a:2U^) 

+ Im(/3)(—£)G(xiu^ + X2U^, a;2U^ — xiu^) 

+ Im(,d)(-£)G(x2U^ - xiu^, + X 2 u‘^) + o(2), 

where 

0(2) = o{xl + xi) + 0(1 Re/3{R)\{x‘l + x|)). 


( 2 . 20 ) 
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Here C = L\e 2 is the projection of the linear operator L onto E 2 . 

Let us write 

J 

Note that for an eigenvector of L* corresponding to by orthogonality relation 

(2.8), we have 

{{-Cf + 4Im(/3)2)(-£)$, u;^) = ^ -4 >j/3r(/ 3| + 4Im(/3)2)(uj, 

K 

= 


where 


rriK = -ISkWk + 4Im(/3)^)(uii-,u|^). 


By (2.14) and the center manifold formula (2.20), we have the following approximation for 
the center manifold 

( 2 . 21 ) 


^ ‘hofcUofc + ^ 4>2m,A:U2m,fc + o(2). 


Here 

= ^l,Kxl + ^2 ,KXiX 2 + ^3,Kxl. 

Using these results, the coefficients of the center manifold function can be computed as 
follows: 


( 2 . 22 ) 


{{Pk + 2Im(/3)^)ff]^^ + 21m{pfgj^ + lm{(3)l3K{gK + 9 k)) , 
$ 2 ,A = {l^xigR + 9k) + 2Im(/3)^ii'(5(f - g]^))) , 

^ 3 ,A = {Wk + 2Im(/3)^)fff + 2lm{j3)‘^g]} - lm{l3)/3K{gK + 9 k)) ■ 
Using ( [2H8l ), ( [2H9l ) and ( [23^ , we find that 

*^1,0^; = ‘hs^ofc) ‘h2,ofc = 0, 


$ 




2, 

= 2 ^ 2 ,2m,k = —^3,2m,k 


Hence the center manifold function (2.21) becomes 

(2.23) ^ = X] ^bOfcUofc(a;? + xj) + 4>i,2m,fc(a;? - 2ixiX2 - xl)u2m,k + o(2). 

k 

Note that 4>i ofc is real while 4>i 2 m,fc is complex. 

Step 3. Construction of adjoint modes. Now we construct the adjoint modes U*^ and U*^ 
orthogonal to and u^. Let us denote the real and imaginary parts of the critical adjoint 
eigenvector by u*^ = ReuJ^^^ ^ and u*^ = ^ and define 


(2.24) 


^,1 ^ (u^u*^)u*^ + (u^u 


U-' 


(u^, U*^)2 


+ (u^u 


1 „*2'\2 


^,2 ^ -(ui,u*2)u*i + (ui,a 


U 


ul u *^)2 + ( u ^, u * 2)2 


Noting that for any two functions of the form fi{x,y) = e' 
a nonzero integer, we have 


_ ^zamTTX, 


9 i{y) for z = 1, 2 where m is 
(Re/i,Re/ 2 ) = (Im/i,Im/ 2 ), (Re/i,Im/ 2 ) = -(Im/i, Re/ 2 ). 
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Thus 

which implies that 
(2.25) 


(u2,U*^) = (u\U*2) = 0, 
(u\U*^) = (u2,U*2) = 1. 


Step 4- Derivation of the reduced equations. Now we write 

(2.26) u(x, y, t) = xi(t)u^(x, y) + X 2 {t)u^{x, y) + $(x, y, t) 
where is the center manifold function, xi,X 2 € M. Note that 

L(u^ + iu^) = /3(u^ + m^) 

implies 

Lu^ = Re(/3)u^ — Im(/3)u^, Lu^ = Re(/?)u^ + Im(/?)u^. 

Also by dehnition of center manifold, we have ($, u*^) = 0 and (<h, u*^) = 0 which by (2.24) 
implies that 

(2.27) (4>,U“) = 0, i = l,2. 


Hence by (2.25) and (2.27) 


(2.28) 


(Lu, U*^) = {xiLu\\J*^) + {x2Lu^, U*^) = Re(/3)xi + Im(/5)x2, 
(Lu, U*2) = (xiW, U*2) + {x 2 Lu^, U*2) = Im(^)xi + Re(/5)x2. 


Plugging ( 2.26[ ) into (2.3), taking inner product with U*® {i = 1,2) and using (2.25) and 

(2.28), we can write the reduced equations as 


= - Im(^)xi + Re(/3)x2 + ((^(u, u), U*^). 
G,(u\u^U*'^) = 0, i,j,A: = l,2, 



dxi 

(2.29) 

dt 

dX2 


dt 

Noting 


and 



G,(4>,$,U*'') = o(3) 


we can expand the nonlinear terms of xi and X 2 in (2.29) 

(2.30) (G(u,u),U*^) = xiG,(ui,4>,U*^') +X2 G,(u 2,$,U*^') +o(3). 


By ( I2.23D , 

(2.31) ^^(u®, 4>, U*^) = ^ 4>i,ofc(a;? + “ 2ixiX2 - 

k 

where we define 

cg'^ = G,(u\uofc,U*^), 

4in,fc = G.(u\u2^,fc,U*^). 


(2.32) 













DYNAMIC TRANSITIONS OF QUASI-GEOSTROPHIC CHANNEL FLOW 


9 


As in (2.17), we can show that 

(2.33) 

Moreover the calculation 

(2.34) 


„ii _ „22 
'-Ok — '-OfcJ 

,11 _ ^22 


^12 _ _„21 
'-Ok — '-Ok-: 


.. ^ 12 ^ 21 
'-2m,k '-2m,k: '"2m,fc '"2m, fc' 


G,(u\uL,U*1) = G,(u\uL,U*2), 

G,(u1,uL,U*1) = -G,(u\uL,U*2) 


implies that 
(2.35) 


'-2m,k '’'-2m,k 


Using (2.30)-(2.35), we have 

(G(u, u), U*-^) = + bf! 2 ixlx 2 + + o(3), j = 1, 2, 

where simple calculations show that 


,2m,kC2m,k: 


(2.36) 


^30 — ^12 — ^21 — ^03 — ^ ^l,OA:Cofc + 

k 

bio = ^12 = -^21 = -^03 = Y1 ^hOfcCofc + i^l,2m,fc4m,fc- 


Thus the reduced equations (2.29) become 

^ = Re(/3)a:i + Im(/ 3 )x 2 + bl^xiixl + xl) - blQX 2 {xj + xl) + o(3), 

(2.37) 

- Im(/3)xi + Re(/ 3 )x 2 + hlQXi{x\ + xl) + blQX 2 ix‘l + xl) + o(3). 
Step 5. Computation of the transition number 7 . Letting z = xi ix 2 , (2.37) becomes 

(2.38) 


f^=/3z + bz\z\^ + oi\z\^), 


where by (2.35) and (2.36) 
(2.39) 


^ — ^30 + *^30 — ^ ‘^1,OA:(CoA: + iCok) 


If Re5 < 0 then (2.38) has a stable limit cycle 


z = re 


iujt 


for Re/3 > 0, i.e. for R > Rq with 


r = \ - 


Re/3 
Re b 


, oj = — Im /3 + Im ~ — Im /3. 


If Reft > 0 then (2.38) has an unstable limit cycle for R < Rq. 

Thus the transition is determined by the sign of real part of ft in (2.39) at R = Rq defined 


as 


(2.40) ^ = Re(ft) = ^ $1,0^4 

since ‘hi ofc, CqI, Cg| are real numbers. 


k=l 
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Step 6. Derivation of the transition number^ in (2.10). 


Using ( |2.18[ ) in ( |2.22[ ), we get 
(2.41) 




all 


-^Ofc(Uofc,U*^) 


The mode uq^ = Ug^ = 0) = {DYQk, 0) corresponds to the eigenfunction with m = 0 in 


(2.6), i.e. solutions of ED^Y = (3D^Y with the boundary conditions T(±l) = D^Y{±1) = 0. 


These solutions are easily obtainable. 

l3ok = - 


k^TT^ 


E, 


(2.42) 


Thus 


yok = 


cos 


sin 


kny 

kny 


if k is odd, 
if k is even. 


2 

(uofc,Uofc) = - / 

a J-i 


DYokV = 


aE 


To obtain 7 , we need to compute 


.11 

-o,fc 


= G,(ui,uofc,U*i) = G(u\uofc,U*^) + G(uofc,ui,U 


*i\ 


(2.43) 
and 

(2.44) gll = G{u^,u^,nl^) 

where is the real part of the critical eigenfunction (with zonal wave number mo) as 


given in ( 2 . 12 ), U*^ is given by (2.24) and G is the trilinear operator (2.13). So we plug 


in 'i/’moi = and into ( 2.43[ ) and (2.44). After tedious 

computations, we can obtain 

aemoTT 


(2.45) 


where 


- - 
%k — 


9ok — 


w 

empirPok 

E 


lm{JiJ2{k)}, 

Im J 3 (/c), 


(2.46) 


Ji = / y*((amo 7 r)^ - D‘^)Ydy, 


-1 


Mk) = / DYokY^ 


( (amovr)^ 


2^2 


k^TT 


- Ydy, 


-1 


Jsik) = / YokYDYdy. 


-1 


Finally we note that Y = Ym^i is either an odd or an even function of y which follows from 
the fact that the equation (2.6) is invariant under the change of variable y —)• —y. Thus 


must be odd (by (2.42), k must be even) otherwise Jsik) is zero. So the nonzero contributions 























DYNAMIC TRANSITIONS OF QUASI-GEOSTROPHIC CHANNEL FLOW 


11 


to 7 in (2.40) comes from even k. Thus we obtain (2.10) from (2.40), (2.41), (2.42), (2.45) 
and (2.46), which concludes the proof. □ 


3. Numerical Evaluation 


3.1. Legendre-Galerkin method for (2.6) and (2.7). In this section we present a method 
to approximate the solutions of the eigenvalue problems (2.6) and d^ . There are fourth- 
order problems so a Legendre-Galerkin method naiis] will be efficient and accurate. 

We look for an approximation of Y in the space = {u G Pn '■ '^(4=1) = u''(±l) = 0}, 
where Pn is the space of polynomials with degree less than or equal to N. 

Using the approach in m , we set 

4 

fjiy) = Lj{y) + '^Cj^kLj+k{y), 

k=l 

with Cj^k to be chosen such that 

/,(± 1 ) = /"(± 1 ) = 0 . 

It is easy to determine from the properties of Legendre polynomials that 

Cj,l = Cj,3 = 0, 

^ 2(2j + 5)(j^ + 5j + 9) 

(j + 3)(j + 4)(2j + 7)’ 

Cj,4 = — 1 — 

and we have Xjv = span{/j : j = 0,1, • • • , N — 4}. 

Writing Y^{y) = Vjfjiy) ^ and plugging it into (2.6), and taking inner product 


with fkiy), {k = 0,1,..., N — 4:) we obtain the Legendre-Galerkin approximation of (2.6) in 
the following matrix form: 


(3.1) 


where 


( R R \ 

E{Ai — 2a^A2 + ~ ~ ^s) j y 

= /3iA2-alAs)y^, 

= {D fj^fk)i ®2j7c = {D fjifk)i ^3,jk — ifjifk)^ 

®4,jfc — (cos Try fj^fk)^ 0'5,jk — {cOSTTy fj, fk) 

Ai — (UjJA:)j,fc=0,...AT—4 ^ — 1, ... 5 


N 


y^ = 


iT 


2/0 


2/1 


2/V-4 


By using the properties of Legendre polynomials, we find that Ai, A 2 , A^ are real symmetric 
banded matrices given by 

r (2 + 2j)(2 + j)(3 + 2j)^(5 + 2y) 

I 0, otherwise 


if j = A; 
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®2 jfc 


- (2 + 2j)(2 + j)(3 + 2j) 

(3 + j)(4 + j) 

< (4(3 + 2j)(5 + 2j)(102 + llOj + A7f + 10 / + j^) 
(3 + j)2(4 + j)2(7 + 2j) 

lo, 


if j = k ±2 


if j = k 
otherwise 


^3,jk 


2(l+j)(2 + j)(3 + 2j) 

(3 + j)(4 + j)(7 + 2j)(9 + 2j)’ 

-8(222 + 196j + 77f + Uf + j^) 

^ (3 + i)(4 + j)(5 + j){6 + j)(ll + 2?) 

4(5580 + 11202j + 9263j2 + 4170/ + 1105j'^ + 168j® + 12f) 

(3 + jT(4 + jf{l + 2j)(7 + 2j){9 + 2j) 

lo, 


if j = A: ± 4 

if j = A: ± 2 

if j = A; 
otherwise 


To approximate A 4 and we fix some integer M and compute the Legendre-Gauss-Lobatto 
quadrature points yn and weights ojn for n = 0,... ,M. Then we compute the matrices 


D‘^fj{yn) and fj{yn)- 


M 

= '^cos{Try)D'^ fj{yn)fk{yn)u}n, 

n =0 

M 

®5jA: — ^ ^ OOs('^?/)/j(?/n)/fc(?/n)^n- 
n=0 


M is chosen large enough to provide sufficient accuracy in the computation of A 4 and A 5 . 
We note that A 4 and A^ are real, full matrices. A^ is symmetric while A 4 is non-symmetric. 


The eigenvalue problem (|3.1|) can be easily solved by using a standard eigenvalue solver. 

Ik if 


)N 


j < k) of (3.1), and corresponding eigenvectors 


For each m G Z, we can numerically find N — 3 eigenvalues (with Ke/3^j > Re/3, 

,j = l,...’N-3. 

By taking the complex conjugate of (3.1), we find = y-m,j^ 
particular (3^^ and 'ipQj are real. 


A 


N 

m,j 


= ^’^d in 


oN 


It is known that (IE]), the computed eigenpairs (/S^^ 
eigenpairs {Pmj, i’mj) of (2.4) exponentially as 00 for 0 < j < 2N/tt for each fixed m 


of ( |3.1[ ) converge to 


Finally, the analog of approximating equation (3.1) for the adjoint problem (2.7) is the 


adjoint problem of (3.1). Namely 


(3.2) 


R 


R 


E(Ai — 20^^42 + a^^s) — /Q:m(^^4 H — ^(vr — a^)A^ — A3) ) y 


O I o' 

TT'^ TT"^ 




= l3*{A2-alA3)y 


*N 


3.2. Numerical computation of the transition number 7 . To approximate 7 in (2.10), 
we follow the following steps: 

Step 1 The parameters of the system are a, e, E and R. We fix the Ekman number E and the 
length scale a. Then by (1.2), the Rossby number e = R/E is also fixed. Moreover, 
to study the transition, R has to be fixed to R = Rq where Rq is determined by 
Assumption 
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E = 0.04 




Figure 1. Neutral stability curves which are defined by Re/3m,i = 0 in the 
a-R plane for zonal wave numbers m = 1, 2, 3, 4 and for E = 0.04, E = 1. 


Step 2 We determine the critical zonal wave number mo and the critical Reynolds number 
i?o ill Assumption]^ InJ^, it was shown that mo = 1 if \/3/4 < a < ao for 
some CIO < Figure show that the difficulties in proving the Assumption 

for ao < a < \/3/2 is purely technical. Extrapolating this result supported by our 
numerical computations of neutral stability curves shown in Figure we claim that 
Assumption can be satisfied by m only if a < \/3/{2m). 

For such m, we set Re (3^ = 0 in (3.1) and solve for the Reynolds number R to find 
Rq which approximates Rq, is the minimum of such Rq^ and the minimizing 
m value is mo- 

Figure shows that Assumption holds for a simple complex pair of eigenvalues 
for almost all a < \/3/2 except for discrete values of a where neutral stability curves 
corresponding to different zonal wave numbers intersect; Assumption 1 is generic. 

computed, the eigensolutions y^o,! = [yj]fjo^ and = [yjjfjo^ of 


nN 

ll'O.m 


Step 3 


With 

(3.1) and (3.2) are also found. 


Step 4 The Legendre-Gauss-Lobatto quadrature points yj and weights LOj are calculated for 
j = 0,..., M where M will be chosen large enough to allow sufficient accuracy in the 
computation of (3.3). 
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Step 5 The values of the eigenmodes and their derivatives at quadrature points yj are com¬ 
puted. 

Af-4 iV-4 

vkfkivj), = Y ykfkivj), 

k=0 k=0 

N-4 N-4 

dyZ,M) = Y y^fkiy,), D^Y^oM) = Y ykfkiyj)- 

k=0 k=0 

Step 6 It is easy to see that if we multiply by a complex number c then 7 is multiplied 

by |cp. To find a unique 7 , we normalize so that maxo<j<M = 1. 

Step 7 We can also normalize so that Ii = 1 in ( |2.9[ ). It is easy to see that this scaling 

has no effect on 7 . 

Step 8 Finally we approximate hik), hik) by 
M 


(3.3) 


h{k) « ^cos(A:7ryj)y^^i(?/j) {{amo-nf - - D'^) Y^^^^{yj)ujj, 

j=0 

M _ 

h{k) ^^s\n{knyj)Y^^^^{yj)DY^^^^{yj)ujj. 
j=0 


Obviously increasing M increases the accuracy of approximation in (3.3). In our 
experiments we found that M = 2N gives enough accuracy. 

3.3. Numerical results. In this section, we present the results of our numerical computa¬ 


tions of 7 for different parameter choices. As discussed in Section 3.2, the only parameters 
that need to be varied are the Ekman number E and the length scale a. 

For showing typical results, we consider a mid-latitude atmospheric jet in a zonal channel 
at a reference latitude Oq = 45°N. The dimensional zonal velocity of the background state 
uq in (1.4) has a maximum U/{tt^E). With a typical zonal velocity of [/ = 15 ms“^, we 


limit our numerical investigations to E values between 0.01 to 0.03 which corresponds to 
maximum zonal jet velocities between 16 ms“^ up to 48 ms“^. For a = 0.2, a typical length 
scale of L = 3000 km yields a channel of length 2L/a = 30000 km in the meridional direction 
that corresponds to about 360° in longitude. We therefore consider a values in the range 
0.1 < a < 0 . 6 . 

For this parameter regime we approximate the transition number 7 following the procedure 
in Section |3.2[ The results we find are presented in Table which suggest that 7 is always 
positive. According to Theorem this means the transition is Type-II at the first critical 
Reynolds number Rq. 

We can also compute the period T = 27r/Im/lmo,i of the solution (2.11) where T has 


been non-dimensionalized by 1/{(3qL). For example, with the above choices of L and U, for 
a = 0.20 and E = 0.01 we find that niQ = 2, Rq = 5.64. The planetary vorticity gradient /3o 
at 45N is 1.6 x 10“^^. We compute the period to be about T = 196 days. The stream function 
of the time periodic solution ( | 2 . 11 [ ) bifurcated on R < i?o which is unstable because of the 
Type-II transition is shown in Figure The pattern indicates a typical one due to barotropic 
instability which, due to the background state zonal velocity, propagates eastwards. 
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Figure 2. The stream function of the bifurcated periodic solution on R < Ro 
at a = 0.2 and E = 0.01; T denotes the period of the periodic orbit. 

Table 1. The value of 7 for 0.1 < a < 0.6 and 0.01 < E < 0.03. 


a 

0.100 0.200 0.300 0.400 0.500 0.600 


0.010 

0.197 

0.197 

0.188 

0.197 

0.194 

0.175 

0.015 

0.174 

0.174 

0.163 

0.174 

0.172 

0.155 

0.020 

0.156 

0.156 

0.143 

0.156 

0.156 

0.143 

0.025 

0.144 

0.144 

0.130 

0.144 

0.146 

0.135 

0.030 

0.135 

0.135 

0.121 

0.135 

0.139 

0.130 


4. Summary and Discussion 

In this paper, we have extended the results of [3] on the existence of a Hopf bifurcation 
in the forced barotropic vorticity equation ( 1 . 1 ) by rigorously showing the type of finite 
amplitude solutions which can occur near this Hopf bifurcation. Central in the analysis is the 
computable quantity 7 which characterizes the transition behavior near the critical point. 

The aim of this paper was only to focus on the theory and numerical computation. As 
an illustration, we explore this number numerically in a parameter regime relevant for an 
atmospheric jet and find that a catastrophic transition is preferred. For typical ocean cases. 
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for example in western boundary currents such as the Gulf Stream and the Kuroshio, and 
the Antarctic Circumpolar Current, the results will be reported elsewhere. 

While (1.1) is a cornerstone dynamical model of the ocean and atmospheric circulation, it 
of course represents only a limited number of processes. As practiced by the earlier workers 
in this field, such as J. Charney and J. von Neumann, and from the lessons learned by the 
failure of Richardson’s pioneering work, one tries to be satisfied with simplified models ap¬ 
proximating the actual motions to a greater or lesser degree instead of attempting to deal 
with the atmosphere/ocean in all its complexity. By starting with models incorporating only 
what are thought to be the most important of atmospheric influences, and by gradually bring¬ 
ing in others, one is able to proceed inductively and thereby to avoid the pitfalls inevitably 
encountered when a great many poorly understood factors are introduced all at once. 

However, the same results will be true if we work on the barotropic equations in primitive 
variables; see (2.2) and [l2]. Second, we would expect that the same type of results obtained 


in this article will also be true if we use higher order friction (e.g. hyper-friction, [13]) in (1.1) 
as this would only change the eigenfunction structure slightly. Third, the method presented in 
this paper combined with the methods introduced in mm can be used to study the case with 
bottom topography |3|, and the case where noise represents ocean eddies. We will explore 
these new directions in the near future. 
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